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Abstract 

Far-field  boundary  conditions  for  external  flow  problems  have  been  developed  based  upon 
long-wave  perturbations  of  linearized  flow  equations  about  a  steady  state  far  field  solu¬ 
tion.  The  boundary  improves  convergence  to  steady  state  in  single-grid  temporal  integration 
schemes  using  both  regular-time-stepping  and  local-time-stepping.  The  far-field  boundary 
may  be  near  the  trailing  edge  of  the  body  which  significantly  reduces  the  number  of  grid 
points,  and  therefore  the  computational  time,  in  the  numerical  calculation.  In  addition  the 
solution  produced  is  smoother  in  the  far-field  than  when  using  extrapolation  conditions. 
The  boundary  condition  maintains  the  convergence  rate  to  steady  state  in  schemes  utilizing 
multigrid  acceleration. 
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1  Introduction 


Numerical  solution  schemes  for  nonlinear  flow  equations  usually  require  the  introduc¬ 
tion  of  one  or  more  artificial  boundaries  to  be  placed  at  some  distance  from  a  body 
around  which  (in  external  flows)  the  flow  takes  place.  On  these  boundaries  appropri¬ 
ate  boundary  conditions  must  be  established  that  are  not  only  physically  correct,  but 
that  also  comply  with  mathematical  requirements  as  well.  In,  for  example,  subsonic 
flow  regimes,  it  is  well  known  that  for  hyperbolic  flow  equations  (e.g.  Euler  Equa¬ 
tions),  one  characteristic  points  back  into  the  flow  domain  at  outflow,  and  therefore 
one  condition  must  be  prescribed;  the  rest,  called  numerical  boundary  conditions, 
must  be  in  some  way  consistent  with  the  partial  differential  equations.  This  paper 
deals  with  the  derivation  of  outflow  (far-field)  boundary  conditions.  These  boundary 
conditions  will  be  both  temporally  and  spatially  local. 

The  proper  formulation  of  far-field  conditions  remains  a  vexing  problem  to  date, 
and  many  authors  have  tried  to  tackle  it  in  a  variety  of  direction's.  It  was  once  common 
in  steady  state  calculations  to  set  p  =  poo,  {p  being  pressure)  and  to,  say,  extrapolate 
all  other  quantities.  Certainly  this  procedure  is  inappropriate  for  boundaries  close 
to  the  body,  and  so  alternate  approaches  are  needed  to  be  developed.  At  present, 
the  Navier-Stokes  solver,  [9],  uses  extrapolation  on  all  physical  quantities  in  subsonic 
regimes  {note  that  in  supersonic  regimes,  extrapolation  of  all  variables  is  correct  from 
a  mathematical  standpoint)  with  great  success.  In  fact,  its  authors  claim  that  they 
have  been  unable  to  replace  far-field  extrapolation  with  another  far-field  condition 
and  to  achieve  satisfactory  steady  state  convergence  [12]. 

It  has  become  popular  to  place  “non-reflecting”  boundary  conditions  at  outflow 
in  subsonic  flows.  While  it  is  true  that  generally  one  cannot  develop  local  boundary 
conditions  that  give  no  reflections,  one  can  consider  conditions  which  to  some  degree 
are  better  than  others.  The  notion  of  better  in  this  case  implies  that  (see  [10]) 

•  the  reflections  are  decreased  and  steady  state  convergence  is  accelerated. 

•  the  reflections  will  decrease  as  the  position  of  the  artificial  boundary  tends  to 
infinity. 

•  as  the  incident  wave  is  more  normal  to  the  artificial  boundary,  the  reflections 
decrease. 

In  1988  Abarbanel,  Bayliss  and  Lustman  (ABL)  [1]  developed  a  non-reflecting 
outflow  boundary  condition  for  viscous,  compressible  flows  over  a  flat  plate.  Their  idea 
was  to  linearize  the  Navier-Stokes  equations  around  a  far-field  steady  state  solution 


Figure  1:  External  flow  topology  for  compressible  viscous  flows  around  NACA012 

with  an  assumed  u-velocity  profile  based  upon  an  asymptotic  approximation  to  the 
equations  themselves.  These  linear  perturbation  equations  were  then  assumed  to  take 
on  a  dispersive  wave  solution  and  the  resulting  set  of  equations  (similar  to  the  Orr- 
Sommerfeld  system)  described  the  asymptotic  behaviour  along  the  outflow  boundary 
-  the  artificial  boundary.  Abarbanel,  Bayliss  and  Lustman  then  assumed  that  the 
decay  rate  of  the  perturbations  is  mainly  controlled  by  the  iteration  scheme’s  inability 
to  dissipate  long  wave  disturbances.  So  they  expanded  the  system,  including  the  wave 
frequency,  around  small  wave  numbers  i.e.  long  waves.  The  resulting  zeroth  and  first 
order  eigensystems  are  solved  for  their  eigenvalues.  The  eigenvalue  of  the  zeroth  order 
eigensystem  physically  corresponds  to  the  decay  rate  of  the  long  wave  disturbances, 
and  the  eigenvalue  of  the  first  order  eigensystem  physically  corresponds  to  the  group 
velocity  of  these  disturbances.  The  slowest  decay  rate  is  found,  and  it  defines  uniquely 
the  group  velocity.  The  details  of  the  aforementioned  strategy  may  be  found  in  [1]. 
The  outflow  boundary  condition  derived  from  this  takes  on  the  form: 

dp  ^dp  ,  . 

^  +  =  (1.1) 

where  a  >  0  is  proportional  to  the  decay  rate  of  the  disturbances,  and  /d  <  0  is  their 
group  velocity.  This  condition  is  very  effective  in  the  flat  plate  case. 

This  paper  uses  the  ABL  approach  in  developing  a  far-field  nonreflecting  boundary 
condition  for  two-dimensional  wake  flows.  We  consider  two  cases.  The  first  being 
compressible  viscous  laminar  flows  around  NACA0012,  as  illustrated  in  figure  1.  In 
this  case  the  compressible  Navier-Stokes  equations  are  numerically  integrated  using 
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Figure  2:  The  flat  plate  topology  for  incompressible  viscous  boundary  layer  flows 

a  time  consistent  scheme.  Multigrid  acceleration  has  been  used  in  this  case  as  well. 

The  second  case  considered  is  incompressible  viscous  boundary  layer  flow  over  a 
finite  flat  plate  as  illustrated  in  figure  2.  In  this  case  the  incompressible  boundary 
layer  equations  are  numerically  integrated  both  time  consistently  and  with  local-time 
stepping. 

The  nonreflecting  boundary  condition  has  been  developed  accordingly  for  time- 
consistent  cases  with  and  without  multigrid  acceleration  and  for  local-time-stepping 
cases  without  multigrid  acceleration. 

In  §2  we  outline  the  long  wave  asymptotic  expansion  procedure  in  the  general 
compressible  case  and  develop  the  appropriate  eigensystem  for  this  case.  In  §3  we 
consider  the  incompressible  boundary  layer  case  and  extend  the  boundary  condition 
to  local-time-stepping  schemes.  Then  in  §4  we  derive  the  non-reflecting  boundary 
condition  both  with  and  without  multigrid  acceleration.  In  §5  we  present  numerical 
results  beginning  with  the  numerical  procedure  used  to  obtain  the  eigenvalues  needed 
in  building  the  boundary  condition,  and  ending  off  with  a  general  discussion  of  the 
results  of  all  the  numerical  tests  performed.  In  §6  we  comment  on  the  results. 
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2  Long  Wave  Asymptotic  Expansions 


The  Navier-Stokes  equations  governing  two-dimensional  compressible  flows  may  be 
written  in  the  following  form: 

continuity  :  pt  -f  -|-  {pv)y  =  0 

^  X  -  momentum  :  p{ut-{- uu^^  +  vUy)  +  p:^  =  ij,(u^^  +  Uyy  +  ^u^c  +  Vy)^)  .  . 

y  -  momentum  :  p{vt  -|-  uv^  +  vvy)  +  Py  p{vxy  -f  Vyy  -|-  +  Vy)y)  ^  ' 

^  energy  :  +  piux  +  Vy)  =  V  ■  (/cVT)  -|-  $ 

where  p  is  the  density,  u  and  v  are  the  x  and  y  components  of  the  velocity  and  p  is 
the  pressure.  In  addition,  /c,  the  heat  conductivity,  Cp  and  c„,  the  specific  heats  at 
constant  pressure  and  volume,  and  their  ratio,  7  =  Cp/c„  are  assumed  to  be  constant. 
The  viscous  dissipation  function,  $  is  defined  as: 

$  =  -f  Vyf  +  2fx{ul  vl)  -t-  pl{ux  +  Vyf  (2.2) 

and  we  shall  use  the  Stokes  relation,  A  -|-  |p  =  0.  The  equation  of  state  for  an  ideal 
gas, 

T  =  (2.3) 

Cp  Oy  P 

is  used.  We  now  define  the  Prandtl  number,  Pr  =  ^,  which  will  be  used  in  what 
follows.  Also  the  Reynolds  number  with  respect  to  the  length  of  the  airfoil  is  defined 
as 

R,  =  (2.4) 

p 

where  the  subscript,  00  implies  that  the  “free-stream”  value  of  the  quantity  is  taken. 
In  addition,  the  viscosity,  p  is  assumed  to  be  constant. 

Assume  that  the  flow  is  external,  and  the  length  of  the  airfoil  is  1.  Defining  the 
dimensionless  variables:  s  =  xf I  and  z  —  y/l,  &n  approximate  wake  profile  in  the  zero 
lift  case  can  be  obtained  [3]: 


^  1  V^i 

=  1  —  rn——p= 

b'oo  2y/'KS 


a 

Ts 


(2.5) 


Integrating  (2.5)  we  obtain  an  expression  for  m  in  terms  of  the  drag  coefficient. 
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We  substitute  the  above  into  (2.5)  and  make  the  following  parabolic  approxima¬ 
tion: 

\  1  if  >  Sf  '  ' 

where  xq  =  sqI  is  far  enough  from  the  airfoil  so  that  the  value  of  m  using  (2.5)  will 
not  change  as  s  >  Sq.  Note  that  (2.7)  satisfies  (2.6),  and  defines  an  approximate  wake 
thickness  of: 

c  3  firs 

Perturbing  (2.1)  in  the  region  <  6^  around  steady  state  conditions  far  down¬ 
stream  gives: 

'  Pt  +  My)Px  +  PooK  +  PooVy  =  0 
PooK  +  MvWX  +  («0)2/U']  +  p;  =  +  U'yy]  +  f  [<  -1-  V'y]^ 

<  PooWt  +  U0(t/)U']  +P'y=  +  V'yy]  +  f  «  +  V'y]y  (2.8) 

p't  +  IPooU'^  +  IPooV'y  +  uo{y)p'x  = 

T^Wxx  +  p'yy)  -  ^{Pxx  +  P'yy)]  +  2  (t  “  l){Uo)ylx[u’y  +  v'^] 

where:  u'  =  u-uq,  p'  =  p-po,  p'  =  p-  po  and  v'  =  v,  and  po  =  Poo,  po  =  Poo,  vo  =  0, 
uo{y)  is  defined  in  (2.7).  In  addition,  ^{•)o  =  ^(Oo  =  0- 

Let 

^  =  (2.9) 


uoiy)  =  U^[K  +  {1  -  K)^] . 


(2.10) 


Note  that  in  the  wake  region,  0  <  A"  <  1  and  outside  of  the  wake  region.  A'  =  1. 
Also,  we  make  the  following  dispersive  wave  ansatz: 


'  P']  [  ^i(j/)  ' 

^  _  .i(;>pt+bx)  ■^2(y) 

F3{y)  • 

L  P'  J  L  F,{y) 

Next,  we  non-dimensionalize  in  the  following  manner: 

Ai  =  Poo(?l  F2  =  U^G2  F3  =  UooGs 

F,  =  p^UlG,  ^  =  ^  b=l 

^  ~  T^s  y  =  ^  =  7^ 


(2.11) 


(2.12) 
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where  u  =  fifpoo  (kinematic  viscosity).  Substituting  the  ansatz  and  the  non-dim- 
ensionalizations  into  (2.8)  leads  to  the  following  set  of  dimensionless  perturbations 
(similar  to  the  Orr-Sommerfeld  system)  in  the  wake  region  (— 1  <  ?/  <  1). 


i[uj  +  /3K(1  +  (^/?2)](?i  +  i^G2  +  =  0 

i[w  +  +  2(1  -  A'),G3  +  i?G,  =  s[G'i  +  if GJ  -  f 

i[ic  +  /3A-(1  +  (1^v^)]G3  +  G;  =  e[lG'i  +  if  Gi  -  /J^Gj) 

i\w  +  +  5^  = 

fJG'l  P\G>  -  Jr)]  +  4(7  -  1)(1  -  K)er,[G',  +  i/SGJ 


(2.13) 


where  (■)'  = 

We  now  perturb  (2.13)  around  long  wavelengths  -  about  small  /?.  To  this  end,  we 
substitute: 

Gi  =  Gf  +  /9Gf >  +  /?%!">  +  •  •  (  (i  =  1, 2, 3. 4)  ,  , 

u;  =  ujq  +  +  /3^a}2  +  ■  *  * 

into  (2.13),  obtaining  the  following  zeroth  and  first  order  problems: 

The  Zeroth  Order  Problem: 


91 =  0 

92  +  fl92  —  2(1  —  K)r}Q.ip  =  0 

'  f V"'  =  0 

s!'  -  4(7  -  1)(1  -  K)P,Mlrig',  +  P,n>p'  -  7£=Aff,n[V.”  +  S&V>]  =  0. 
where  gi  =  Gf\  i  =  1,2, 3, 4  together  with  the  change  of  variables: 

iwo  =  — efl  53  =  gt  =  . 


(2.15) 


(2.16) 


The  First  Order  Problem: 

Gi  -  $'  =  ji92  +  ^[a;i  +  /i;'(l  +  ^^ri^)]9i 

G'i  +  Q.G2  -  2^(1  -  K)ri^  =  i[u;i  +  K{1  +  ^r,‘')]92  +  -  y] 

-  y  =  i|[cc;a  +  K{1  + 

G'i  -  4(7  -  1)(1  -  K)PrMlrjG'2  +  = 

-ie^PrMlnico^  +  K{1  + 

I  -iPr92  +  4£S(7  -  1)(1  -  K)PrM^gnip 

with  the  following  change  of  variables: 


(2.17) 


G,  =  £GS'>  Gi  =  £Gf’> 


(2.18) 
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In  order  to  solve  the  zeroth  and  first  order  problems,  appropriate  boundary  con¬ 
ditions  need  to  be  established.  We  shall  derive  these  conditions  by  looking  at  the 
perturbations  outside  the  wake  region,  if  =  1,  and  then  match  at  =  ±1.  One 
can  easily  check  that  placing  K  =  I  \n  (2.13)  gives  a  system  with  constant  coef¬ 
ficients.  We  shall  call  the  perturbations  outside  the  wake  region,  Q,  analogous  to 
G  in  the  wake  region.  Clearly  this  implies  that  Qi  takes  the  form  with  Qi 

as  constant.  Now  the  exact  boundary  conditions  for  the  zeroth  and  first  order  sys¬ 
tems  are  Gi(±l)  =  =  1,2, 3,4.  In  order  to  get  a  relationship  on  the  Gj’s 

we  use  the  fact  that  they  satisfy  Since  also  Qi  must  satisfy 

^  i  =  1,2, 3,4,  then  clearly  appropriate  boundary  conditions  on 

Gi(±l)  are: 

^  i  =  l,2,3,4.  (2.19) 

/  J  r)=±l 

In  order  to  obtain  r  in  the  previous  expression,  we  do  the  following: 

Substituting  K  =  1  into  (2.13)  and  using  Qi  =  Qic'^'^  and  hence  Q'^  =  rQi,  and 
Q'l  —  'f'^Qi  we  obtain: 


i(u)  -f-  P)Qi  -|-  i^Q2  +  fQz  =  0 
{i{u  +  15)-  er^)Q2  -  i^^rQ^  -1-  i^Q^^  =  0 
(i(a;  -f-  ^)  -  4Ser^)Q3  -  r(52  +  rQ4  =  0 

(i(a;  -t-  ^)  -  ^^er^)Q4  +  p;M^r^Q2  +  j^(?2  + 


(2.20) 


-Q3  =  0. 


We  expand  Qi  and  r  around  ^  as  in  (2.14)  and  define  =  qi  and  obtain  the 
zeroth  order  problem  in  the  free  stream  region: 


iLOoqi  +  roq3  =  0 

{iujo  -  £rl)q2  =  0 
(iuo  -  ^erDqs  +  roq4  =  0 
+  ^093  +  (iwo  -  rlj 


(2.21) 


-+)Mlq4  =  0 


where  it  is  evident  that  either  92  =  0  or  =  ^  =  —  fl.  We  note  that  e’"’^  cannot 
decay  as  t]  ±00  if  Tq  =  —Cl  hence  92  =  0.  For  the  rest  of  the  system,  the  first, 
third  and  fourth  equations  in  (2.21),  we  apply  the  existence  criterion  for  nontrivial 
solutions  of  homogeneous  systems,  namely: 


iu>o  -  |£r^ 


MUicao  -  ^erl) 


(2.22) 
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We  expand  this  around  the  first  row  retaining  terms  up  to  0{e^M^).  Assuming  that 


<C  1  ,  then,  we  find  that; 


1  +  (| 


(2.23) 


For  subsonic  flows  with  7  =  1.4  and  Pr  =  0.72  and  small  e  (physically  reasonable) 
the  latter  result  is  negative.  Hence  the  perturbations  will  not  decay  in  7/.  Therefore, 
we  choose; 

ro  =  q=^^Moo£  (2.24) 

where  the  —  sign  is  taken  for  positive  rj  and  the  +  sign  is  taken  for  negative  7,  so 
that  the  perturbations  in  free-stream  will  die  out  in  77. 

The  first  order  problem  in  the  free  stream  region  takes  the  form: 


=  -*(wi  +  1)91  -  riqs 
{iojo  -  erl)Q^2^  =  i^roqz  -  iq^ 

{iujo  -  +  1)  +  feJ'orOga  -  ng4 

+  roQ^z^  +  {wo  -  = 

(-7(071  +  1)  +  ^£rori)M^94  -  ^^qi  -  nqz  ■ 


(2.25) 


The  To  chosen  above  makes  the  coefficient  determinant  of  the  left-hand-side  singular. 
Hence,  solvability  here  requires  that  the  same  coefficient  determinant  with  one  of  its 
columns  replaced  by  the  right-hand-side  be  singular  as  well.  This  criterion  leads  to; 


ri  =  ±7(a;i  -f  l)Moo  +  0{e‘^Ml,). 


(2.26) 


In  order  to  fix  the  boundary  conditions  for  (2.15)  and  (2.17),  we  recall  the  require¬ 
ment  that: 

-j—  =  1)2, 3,4 

/  -I  77=±1 

whose  ^  expansion  is 

1  [g!»'  +  /3G!-'  +  . .  1  =  [(r„  +  l}r,  +  ..  .)(gS»'  +  ^Gf>  +  . . 


gives 

f  G!'’>'(±1)  =  r„GS“> 

I  G!’>'(±l)  =  r„G!"(±l)  +  r,Gr>.  ' 

This,  together  with  $2  =  0  •«->■  ^2(^1!)  =  0  provides  us  with  appropriate  boundary 
conditions  for  the  zeroth  and  first  order  perturbation  problems.  The  reader  will  note 
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that  the  systems  (2.15)  and  (2.17)  are  of  order  six,  so  only  six  conditions  will  be 
prescribed  for  them.  For  physical  reasons,  we  do  not  place  conditions  on  the  pressure 
perturbations,  i.e  on  xf)  and  on  II. 

In  summary,  the  boundary  conditions  on  (2.15)  are: 

^((±1)  =  ro5i(±l) 

g2{±l)  =  0  (2.28) 

^'(±1)  =  ro¥?(±l) 

and  on  (2.17)  are: 


=  roGi(±l)  -t-  £ri5fi(±l) 

^2(^1)  =  ro(72(±l)  (2.29) 

$'(±1)  =  ro$(±l)  +  £rii^(±l) 

Since  the  perturbations  decay  like  we  solve  the  zeroth  order  problem  for  its 

minimal  positive  fl  and  then  use  this  U  in  the  first  order  problem  to  find  wi  using  a 
Fredholm  like  process.  Since  ^  l/3=o=  5  then  ha.s  the  physical  meaning  of  being 

the  group  velocity  of  the  longest  waves. 


3  The  Incompressible  Case 

In  this  section  we  shall  find  appropriate  and  coi  in  the  special  case  of  incompressible 
viscous  flows,  i.e.  =  0. 


3.1  The  Time  Consistent  Case 

The  zeroth  order  problem  in  this  case  is: 


=  0  /  /  I  n  _  n 

g!^  +  ilg,-2{l-K)gu^  =  0  J  , 


g'y  -f  CtPr<p'  —  0 


whose  normalized  solution,  in  the  sense  that  ||  g2  Wjjx  ^=1  is  as  follows: 


=  yj  =  ^0  =  0  g2{T})  -  cos{-g) 


hence: 


0  ^  ^min.  — 
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The  first  order  problem  is  interesting  in  that  looking  at  (2.29)  one  would  think 
that  the  appropriate  boundary  conditions  on  are  ^  making  the  usage 

of  the  Fredholm  alternative  incorrect.  It  can  be  easily  shown  that  in  fact: 


r  Gi  -  $'  = 

G'i  +  ^G2-29.{1-  K)r}d> 
'  $"  +  -  I n' = 

i  G'i^  Prm' = -iPrg2 


+  K  {I  +  ^^''f)]g2 


G[{±1)  =  0 
^2(^1)  —  0 
$'(±1)  =  0 


(3.4) 


so  that  left-hand-null- vector  of  the  first  order  problem  is  indeed  cos[y).  Taking  the 
scaler  product  of  co$(y)  and  the  second  equation  in  (3.4)  gives 

u.,  =  -K-{l-  ^)(i  -  |) .  (3.5) 

It  should  be  made  clear  at  this  point  that  had  the  condition  5^2 (il)  =  0  been  used 
instead  of  5'2(±1)  =  0,  in  the  zeroth  order  problem,  we  would  have  found  that  g2{v)  = 
sin{y).  In  this  case  we  would  not  have  been  able  to  solve  the  first  order  system  for 
u>i.  So  it  is  very  important  to  identify  the  proper  boundary  conditions  in  order  to 
assure  finding  fi  and  a>i. 

The  longwave  expansion  asymptotics  done  thus  far  are  based  on  the  assumption 
that  the  system  of  partial  differential  equations  are  solved  in  a  time  consistent  manner 
using  a  global  time  marching  scheme.  In  addition,  had  we  started  with  the  viscous 
incompressible  boundary  layer  equations 


Ux  -|-  Vy  -  0 

Ut  -I-  {u^)x  +  {uv)y  =  VUyy 


(3.6) 


and  done  the  same  type  of  longwave  asymptotics  as  in  §2  (note  that  the  far-field 
profile,  (2.7)  is  appropriate  for  this  case  as  well)  we  would  have  obtained  the  same 
fl  and  LO\  as  above.  This  means  that  essentially  we  have  found  the  eigenvalues  for 
the  viscous  incompressible  boundary  layer  equations  for  time  consistent  numerical 
schemes.  In  the  next  subsection,  we  shall  find  these  eigenvalues  when  local  time 
stepping  is  to  be  used. 


3.2  The  Local  Time  Stepping  Case 

Many  researchers  use  local  time  stepping  in  numerical  marching  schemes  towards 
steady  state.  While  local  time  stepping  is  not  time  consistent,  it  has  been  found  to 
accelerate  convergence  to  steady  state. 

We  return  to  the  viscous  incompressible  flat-plate  boundary  layer  equations,  (3.6). 
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We  are  interested  in  integrating  this  system  numerically.  In  order  to  ensure  nu¬ 
merical  stability,  the  Courant  Fredrichs  Levy  (CFL)  condition  on  the  time  step  for 
each  grid  node  requires  (see  [7])  that  At  be  bounded  by  some  function  of  Axjj  and 
Ayij.  The  function  depends  upon  the  numerical  scheme  used.  In  our  numerical  ex¬ 
periments  Runge  Kutta  schemes  were  used.  We  were  conservative  in  our  approach 
and  required  that 

At  <  minA%  (3.7) 


where 


Adij  <  min 


ttij  1 


where  0  <  ^  <  1.  The  addition  of  the  factor  ^  helped  stabilize  the  numerical  iterative 
process.  For  the  runs  with  global  time  stepping  on  a  grid  with  constant  Ax  and  Ay, 
0  =  1  was  enough  for  stability,  however  for  the  runs  on  a  grid  with  nonconstant  Ay, 
we  chose  /?  =  |  to  achieve  stability  (see  §5). 

Local  time  stepping  means  that  at  the  node  {xij,  yij)  a  time  step  of  A%  is  chosen. 
Using  local  time  stepping  one  discretizes  the  temporal  derivatives  in  the  following 


manner: 


Suppose  we  discretize  (3.6)  on  a  non-uniform  rectangular  grid  and  temporally 
discretize  as  in  (3.9).  Provided  that  the  spatial  derivatives  are  discretized  consistently, 

signifying  that  the  grid  is  not  changing  in  time.  Therefore, 
these  spatial  discretizations  are  in  fact  valid  approximations  for  the  spatial  derivatives. 
However,  the  temporal  derivative  is  quite  different.  Since 


_ 


we  clearly  see  that  what  we  are  approximating  is  in  fact 


At-*o  \A$ij  J  dt 

Suppose  that  the  grid  was  constructed  so  that  in  the  boundary  layer, 

and  that  the  grid  is  geometrically  expanding  in  the  y  direction  so  that 

Ayij  =  k  >  1 


(3.10) 


(3.11) 
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then 


2 


and  therefore 


lin>  _ 


2  F(V,  k). 


(3.12) 


Note  that  0  <  F{y;  k)  <  1  when  A;  >  1  and  that  F{y,  1)  =  1. 

Therefore,  while  we  are  solving  (3.6)  on  this  special  grid  using  local  time  stepping, 
we  are  in  actuality  solving  the  modified  system: 


Ux  +  Vy  =  0 

F{y,  k)ut  +  {u^)x  +  {vu)y  —  VUyy 


(3.13) 


We  shall  now  carry  out  a  long- wave  asymptotic  analysis  for  the  above  system  of 
equations.  We  note  that  the  far-field  wake  profile,  (2.5)  is  still  appropriate  in  this 
case.  Hence,  using  the  parabolic  far-field  approximation  for  the  wake  profile  (2.7)  we 
again  define  the  primed  variables 


u  =  u  —  Uo  V  =  V  . 


Perturbing  (3.13)  around  steady  state  variables  retaining  linear  terms  in  primed  vari¬ 
ables  and  substituting  the  dispersive  wave  ansatz: 


Fi(y) 

Hv) 


(3.14) 


and  nondimensionalizing  0,  6,  e,  and  y  as  in  (2.12)  we  obtain  the  following  system  of 
ordinary  differential  equations  in  y. 

iC?i(r/)/?-fG''(7?)  =  0 

i[F(rr,  k)L>  +  /JA-(1  +  1^,2)JG,(,)  (3.15) 

+2G2(,)(1  -  if),  =  eG'M 

where  Fi  =  U<x,G\  and  F2  =  U00G2  and 


nv,k)  =  \ 

1 

Perturbing  Gi  and  uj  about  P: 

UJ  =  UJq  -f-  ^UJi  -f  jFuJ2 


\v\  <  1 
\v\  >  1 


(*  =  1,2, 3, 4) 


(3.16) 


(3.17) 
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and  substituting  these  perturbations  into  (3.15)  we  get  the  following  zeroth  and  first 
order  problems  in 


The  Zeroth  Order  Problem; 


<^'(7?)  =^0 

g'(  +  k)gi  -  2(1  -  K)T}Q,ip  =  0 


where  =  gi,  g2  =  eCl(p  and  wq  =  —Q,e. 

The  First  Order  Problem: 


r=-4s- 


Gf  +  fiF(i);  k)Gi  -  2n*(l  -  K)r,  = 
Vi  [rfe;*;)"! +  *■(!  + 


(3.18) 


(3.19) 


where  Gi  =  and  ^ 

Again,  we  need  to  match  solutions  at  \rj\  =  1  (at  the  boundary  of  the  wake 
region)  —  when  K  =  1.  Formally  setting  k  =  1  because  out  of  the  wake  region 
the  grid  stretching  should  have  no  effect  either  on  or  on  cui  we  obtain  the  linear 
perturbation  system  in  the  free  stream  region: 


f  iGi{T])0  +  G2iTi)  =  0 

(  i[u:  +  ^]Gi{rj)  =eG1{ri) 

whose  constant  coefficients  imply  a  solution  of  the  form: 


Giiv) 

G2{v) 


(3.20) 


(3.21) 


Expanding  as  before,  around  ^ 

Gi  =  gi-V^G]^--- 

r  =  To  +  )3ri  d - 

w  =  ujq  “b 

we  obtain  the  following  zeroth  and  first  order  problems  in  the  free  stream  region: 


0*^  order  : 


roh  =  0 

{iujQ  -  erDg-i  =  0 


(3.22) 
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from  which  immediately  ^2  =  0.  Since  iuo  =  then  =  0  as  well,  otherwise  the 
perturbations  will  not  decay  in  tj. 


1^*  order  (I3^)  :  igi  +  rQG\  +  rig2  —  0 

G\{iu}o  -  er^)  =  0 

As  in  the  zeroth  order  problem  above,  both  GJ  =  0  and  G\  =  0. 

This  implies  that  the  boundary  conditions  for  the  zeroth  order  problem  are 


(3.23) 


5'i(±l)  =  0 

<^(±1)  =  0 


(3.24) 


and  that  the  boundary  conditions  for  the  first  order  problem  are 


Gi(±l)  =  0 

$(±1)  =  0  . 

Notice  that  the  zeroth  order  problem  reduces  to: 

g"  +  0,F{t);  k)gi  =  0 
^i(±l)  =  0 


(3.25) 


(3.26) 


Since  0  <  F{r);  k)  <  I  when  A;  >  1  then  the  Sturm  Comparison  and  Oscillation 
theorems  (see  [4])  indicate  that  the  minimal  positive  eigenvalue  17  in  (3.26)  will  be 
greater  than  the  minimal  positive  eigenvalue  in  the  non-stretched  incompressible  case. 
This  fact  is  crucial  and  to  a  certain  extent  indicates  why  local  time  stepping  is  so 
effective.  Since  it  is  believed  that  indeed  the  long  wave  perturbations  are  the  slowest  to 
converge  to  steady  state,  then  we  have  just  shown  that  the  decay  of  these  disturbances 
is  faster  in  local  time  stepping  regimes  than  in  global  time  stepping  regimes. 

In  (§5)  we  present  numerical  evidence  indicating  that  the  nonreflecting  boundary 
condition  is  quite  effective  in  local  time  stepping  regimes  indicating  that  the  absorbed 
long  waves  do  in  fact  decay  as  we  predict.  This  is  an  interesting  illustration  of  one  of 
the  reasons  why  local  time  stepping  is  effective. 


4  Derivation  of  the  Far-Field  Non-Reflecting  Bound¬ 
ary  Condition 

Once  the  two  eigenvalues,  17  and  coi  are  found,  an  appropriate  far-field  boundary 
condition  can  be  developed  to  accommodate  the  physical  situation  at  the  far-field. 
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4.1  Derivation  for  Single-Grid  Numerical  Methods 

From  the  ansatz  (2.11),  we  have: 

Substituting  the  dimensionless  temporal  and  spatial  variables:  into  the  profile,  we 
obtain  (up  to  order  /3): 

P  =  +  Poo  . 

Taking  a  derivative  with  respect  to  0  and  using  the  identity  isojuo  =  we  arrive 

at  the  far-held  boundary  condition: 

dp  4  -  ,  .  dp 

_  =  +  (4.1) 

In  the  incompressible  case  where  fl  =  condition  takes  the  form: 

dp  IT ..  dp  ,  ^ 

^  =  -g(p-p^)  +  wi^.  (4.2) 


4.2  Derivation  for  Numerical  Methods  Utilizing  Multigrid 
Acceleration 


The  non-rehecting  boundary  condition  just  developed  was  based  upon  an  analysis 
around  long  waves.  Should  multigrid  acceleration  be  used  in  conjunction  with  the 
numerical  scheme  then  short  waves  must  be  taken  into  account  as  well.  Substituting 
P  for  p  —  Poo,  the  boundary  condition  (4.1)  may  be  rewritten  as 


dP 

dt 


_4_ 


flP  “t"  Wj 


dx 


(4.3) 


Taking  the  Fourier  transform  in  x  of  the  above  equation  we  obtain 

Pt  =  -^flP  -1-  iku-iP.  (4.4) 

which  has  the  general  solution: 


For  short  waves,  or  high  Fourier  modes,  k,  the  highly  oscillatory  behavour  of  the 
symbol  P  with  t,  will  cause  short  wave  disturbances  in  the  numerical  solution.  This 
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is  especially  critical  when  using  multigrid  techniques  whose  goal  on  each  grid  is  to 
damp  short  waves.  Our  original  nonreflecting  boundary  condition,  while  effective  in 
damping  long  wave  disturbances  caused  by  reflections,  may  create  new  short  wave 
disturbances  and  destroy  the  effects  of  the  multigrid  acceleration. 

By  replacing  ik  in  (4.4)  with 

ik 

1  +  ik 

we  obtain  a  new  Fourier  symbol  free  of  oscillatory  behaviour  even  as  k  gets  large: 

P=  (4.5) 

The  symbol  (4.5)  comes  from  the  equation: 

P,  =  ~aP  +  -pT^P . 

yTT  1  +  tfc 

Transforming  back  to  the  original  variables  we  obtain  the  boundary  condition  to 
be  used  with  multigrid  acceleration: 

dp  4n  dp  d  \dp  4:  1 

A  =—(?-!’«•) +  (4.6) 


d\dp 


dx  I  dt 


^  =  -7ir(p-Poo)+  wi-7rr 


4fi]  dp 


97rJ  dx 


5  Numerical  Results 

In  this  section  we  present  various  numerical  results  obtained  using  the  far-field  bound¬ 
ary  condition  just  derived.  We  break  this  section  into  three  subsections.  The  first 
describes  the  numerical  technique  for  calculating  Vl  and  wi  in  compressible  cases.  The 
second  section  deals  with  the  numerical  solution  of  the  incompressible  boundary  layer 
equations  both  with  global  and  local  time  stepping.  The  third  presents  the  results 
obtained  in  compressible  viscous  flows  around  the  NACA0012  airfoil  using  global  time 
stepping  regimes  with  and  without  multigrid  acceleration. 


5.1  The  numerical  calculation  of  Q  and  uji  in  the  global  time 
stepping  case 

The  discretized  form  of  the  zeroth  order  problem  (2.15,  2.28)  is  a  generalized  eigen¬ 
value  problem: 

Ax  =  h^Q,B{Cl)x  (b.l) 
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where  h  is  the  cell  size,  A??,  and  x  is  {g'^,g^,(p,'tp).  The  problem  is  non-linear  since 
B  =  -B(fi)  is  a  function  of  fl,  so  the  following  iterative  scheme  was  used: 

\  Ax  =  h^a^+^B{n^)x 

where  at  each  iteration,  the  smallest  positive  0,  was  determined.  Generally  conver¬ 
gence  was  obtained  after  three  iterations  and  the  main  result  is  (Figure  3): 

For  all  physically  relevant  values  of  e,  Moo  and  K  (see  eg.  2.9); 


The  discrete  first  order  problem  (2.17,  2.29)  can  be  written  in  the  form: 

C{n)x  =  {A-  h^nB{^))x  =  n^+  (5.3) 

where  X  =  (Gi,  (?2,  $, 'F)  and  where  A  ^  A  and  B  <r^  B  except  that  the  boundary 
conditions  on  g2  are  not  parallel  to  those  on  G2.  The  technique  used  to  obtain 
was  to  find  the  left  hand  null  vector  of  (7(17),  u,  and  then 


(n,:^2) 


(5.4) 


Unlike  17,  wi  was  not  found  to  be  constant  for  all  physically  relevant  parameters.  Fig¬ 
ure  4  graphically  shows  the  values  of  loi  as  a  function  of  1  -  A'  for  Moo  =  0, 0.1, 0.2, 0.3 
(for  neatness  Moo  =  0.4, 0.5  were  not  shown).  The  main  result  is: 
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wi  — >  —1  for  0  <  M  <  0.5  and  for  all  physical  e. 

The  limit  M^o  0  is  singular. 

In  Table  1  values  of  Q,  and  for  a  variety  of  parameters,  s,  Moo  and  K  for  time  con¬ 
sistent  numerical  integration  schemes  are  listed.  These  values  were  obtained  through 
a  program  written  in  MATLAB. 

5.2  Incompressible  Viscous  Flow 

5.2.1  Case  1-  Global  Time  Stepping  Regimes 

The  viscous  incompressible  boundary  layer  equations  (for  a  flat  plate)  (3.6)  have 
been  solved  in  the  geometry  of  Figure  2.  The  Reynolds  number  (with  respect  to 
the  plate  length),  Ri,  taken  was  100000,  and  a  rectangular  grid  (476x61)  with  cell 
aspect  ratio,  of  25  was  used.  The  plate  had  76  nodes  in  the  x  direction,  and 
there  were  400  x  nodes  in  the  wake  region.  At  inflow,  one  quarter  the  way  down  the 
flat  plate,  a  Blasius  profile  was  placed.  The  outflow  boundary  was  placed  four  plate- 
lengths  from  the  trailing  edge.  The  momentum  equation  was  integrated  in  time  using 
a  Runga-Kutta  scheme  as  suggested  by  [6],  where  the  first  order  spatial  derivatives 
were  discretized  using  a  compact  fourth-order  scheme  developed  by  Abarbanel  et. 
al.  [2].  The  second  derivatives  were  discretized  using  standard  second-order  central 
differencing. 
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K 

e 

Moo 

n 

0.950000 

0.000100 

0.100000 

2.462296 

-0.979065 

0.950000 

0.001000 

0.100000 

2.462296 

-0.977492 

0.950000 

0.000100 

0.200000 

2.462190 

-0.979052 

0.950000 

0.001000 

-0.975585 

0.950000 

0.000100 

0.300000 

2.462014 

-0.978915 

0.950000 

0.001000 

0.300000 

2.461997 

-0.973317 

0.950000 

0.000100 

0.400000 

2.461773 

-0.978703 

0.950000 

0.001000 

0.400000 

2.461732 

-0.982999 

0.950000 

0.000100 

-0.977015 

0.950000 

0.001000 

0.500000 

2.461392 

-0.967744 

0.800000 

0.000100 

0.100000 

2.461776 

-0.911913 

0.800000 

0.001000 

0.100000 

2.461767 

-0.906617 

0.800000 

0.000100 

0.200000 

2.460254 

-0.904457 

0.800000 

0.001000 

0.200000 

2.460203 

-0.894001 

0.800000 

0.300000 

2.458007 

-0.897642 

0.800000 

■imiiHiiiiM 

0.300000 

2.457900 

-0.883191 

0.800000 

0.000100 

0.400000 

2.455157 

-0.892715 

0.800000 

0.001000 

msm 

-0.875650 

0.800000 

0.000100 

0.500000 

2.451718 

-0.892979 

0.800000 

0.001000 

0.500000 

2.451536 

-0.870731 

0.700000 

0.000100 

0.100000 

2.461120 

-0.862841 

0.700000 

0.001000 

0.100000 

2.461103 

-0.855154 

0.700000 

0.000100 

0.200000 

2.458011 

-0.847079 

0.700000 

0.001000 

0.200000 

2.457942 

-0.833530 

0.700000 

0.000100 

0.300000 

2.453518 

-0.840599 

0.700000 

0.001000 

0.300000 

2.453408 

-0.820463 

0.700000 

0.000100 

0.400000 

2.447673 

-0.835958 

0.700000 

0.001000 

0.400000 

2.447541 

-0.814361 

0.700000 

0.000100 

0.500000 

2.440383 

-0.834136 

0.700000 

0.001000 

0.500000 

2.440236 

-0.812478 

Table  1:  fl  and  for  various  ii',  e,  and 
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Figure  5:  Convergence  behaviours  for  inconapressible  case. 


Two  cases  were  run.  The  first  with  extrapolation  of  u  at  outflow,  and  the  second 
using  our  far-field  boundary  condition  (incompressible,  (4.2))  at  outflow  with  up¬ 
dated  every  iteration  as  in  (3.5).  In  each  case  5000  time  steps  were  taken.  Figure  5 
shows  the  convergence  behaviour.  It  is  evident  that  the  non-reflecting  condition  re¬ 
duced  markedly  the  residuals  in  the  flow  field,  implying  that  the  the  non-refleflecting 
boundary  condition  produced  a  more  stable  numerical  solution  than  did  extrapola¬ 
tion. 

In  order  to  test  the  effect  of  the  nonreflecting  boundary  condition  as  the  outflow 
boundary  is  moved  closer  to  the  trailing  edge,  additional  runs  were  performed  with  the 
outflow  boundary  at  s  =  2  and  at  s  =  3.  Figure  7  shows  that  as  s  decreases,  the  effects 
of  the  nonreflecting  boundary  condition  are  more  dramatic.  In  addition,  Call  S2,  S3, 
and  S4  the  solutions  of  the  flatplate  problem  on  grids  G2  C  G3  C  GA  respectively.  Use 
N  as  an  extention  for  a  case  with  the  nonreflecting  boundary  conditions  at  outflow, 
and  use  E  as  an  extention  for  a  case  with  extrapolation  conditions  at  outflow.  We 
looked  to  see  to  what  degree 


SN2  c  SNZ  c  5iV4, 


and 


SE2  C  SE3  C  SEL 


It  turns  out  that  SN2  C  SNZ  C  5'iV4  up  to  four  significant  digits  whereas  SE2  C 
SEZ  C  SE4  only  up  to  two  significant  digits.  Hence,  we  are  certain  that  our  bound¬ 
ary  condition  maintains  accuracy  better  than  the  extrapolation  condition  does.  All  of 
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Figure  6:  Residuals  after  5000  time  steps  -  Extrapolation  (top)  Non-reflecting  (bot¬ 
tom) 

this  indicates  that  one  could  use  a  smaller  computational  domain  while  using  the  non¬ 
reflecting  boundary  condition  and  maintain  “large  computational  domain”  accuracy, 
thereby  dramatically  reducing  the  amount  of  computational  work  needed. 


5.2.2  Case  2-  Local  Time  Stepping  Regimes 

We  shall  now  describe  the  results  obtained  using  the  nonreflecting  boundary  condi¬ 
tion  in  conjunction  with  a  local  time  stepping  regime.  Let  us  state  the  conditions  of 
the  case  studied. 

(1)  The  Grid:  As  in  subsubsection  5.2.1,  the  global  time  stepping  case,  the  equa¬ 

tions  are  discretized  on  a  476  x  61  rectangular  grid.  The  inflow  boundary  is  located 
one  quarter  away  down  the  flat  plate  and  the  outflow  boundary  is  placed  4  chord 
lengths  down  stream  from  the  trailing  edge.  While  Aa;  is  constant.  Ay  grows  geomet¬ 
rically  so  that  Aj/ij  =  At/j,i,  where  Ayij  =  The  expansion  factor,  k, 

was  chosen  so  that  in  60  intervals,  y  reaches  one  chord  length.  It  should  be  noted  that 
the  aspect  ratio  of  cells  above  the  plate  and  the  centerline  are  set  identical  to  those 
in  the  global  time  step  case  so  that  in  the  local  time  stepping  case  At  is  determined 
by  the  parabolic  part  of  the  system  in  most  of  the  wake  region,  recall  (3.11). 

(2)  Physical  conditions:  In  this  case,  Ri  is  again  taken  to  be  100000.  Also,  the 
boundary  condition  at  inflow  is  the  Blasius  profile  one- quarter  way  down  the  plate. 
The  outflow  boundary  is  located  4  chord-lengths  away  from  the  trailing  edge.  Equa- 


tions  (3.18)  and  (3.19)  with  their  boundary  conditions,  (3.24)  and  (3.25)  are  numer¬ 
ically  solved  in  this  case  for  Ct  and  wi  using  a  program  written  in  Mathematica.  The 
algorithm  is  similar  to  the  the  one  presented  in  subsection  5.1.  In  this  case,  is  not 
a  constant  and  depends  upon  the  expansion  factor,  k.  For  these  physical  conditions 
we  obtain 


n  =  5.95497 

2 

which  is  significantly  larger  than  obtained  in  the  global  time  stepping  case.  This 
implies  that  the  long  wave  disturbances  decay  faster  when  local  time  steps  are  used 
than  with  global  time  steps.  Regarding  cji,  recall  that  in  the  global  time  stepping 
scheme  case  the  relationship  between  between  uji  and  K  was  linear.  Therefore,  the 
relationship  between  loi  and  cj.  was  linear  as  well.  It  has  been  numerically  determined 
for  the  present  local  time  stepping  case,  that  the  relationship  between  Wj  and  Cd  is 
linear  as  well,  more  specifically: 
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Figure  8:  Residual  plots  after  5000  iterations  (a)  extrapolation,  gts  (b)  non-reflecting 
condition,  gts  (c)  extrapolation.  Its  (d)  non-reflecting.  Its 


wi  =  56.32crf  -  2.01413  .  (5.5) 

The  non-reflecting  boundary  condition  with  these  coefficients  was  used  as  an  out¬ 
flow  boundary  condition,  as  was  extrapolation.  The  results  were  compared  and  we 
state  them  now. 

In  figure  (8)  the  pointwise  residuals  after  5000  iterations  for  the  four  test  cases  (ex¬ 
trapolation  outflow  conditions-global  time  stepping,  non- reflecting  outflow  conditions- 
global  time  stepping,  extrapolation  outflow  conditions-local  time  stepping,  and  non¬ 
reflecting  outflow  conditions-local  time  stepping),  are  plotted.  In  figure  9  the 
norm  of  the  residuals  (for  each  case)  is  plotted  as  a  function  of  time  steps.  Figure  9 
(a)  graphs  the  logarithm  of  the  norm  of  the  residual  as  a  function  of  time  steps 
(iterations)  for  global  time  stepping  runs.  Figure  9  (b)  graphs  the  logarithm  of  the 

norm  of  the  residuals  as  a  function  of  time  steps  for  local  time  steping  runs.  The 
bottom  dotted  line  is  a  special  case  that  will  be  discussed  later  on  in  this  section.  It 
is  clear  from  the  results  that 

•  Local  time  stepping  regimes  converge  faster  than  their  global  time  stepping 
counterparts. 

•  The  case  with  non-reflecting  outflow  conditions  with  global  time  stepping  reaches 
steady  state  faster  than  the  case  with  extrapolation  outflow  conditions  and  local 
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(a) 


(b) 


time  stepping. 

Note  that  adding  the  non-reflecting  boundary  condition  to  the  local  time  stepping 
scheme  speeds  up  the  convergence.  This  implies  that  the  non-reflecting  boundary 
condition  accelerates  convergence  to  steady  state  above  and  beyond  local  time  step¬ 
ping. 

For  experimental  purposes  only  we  ran  the  flat  plate  code  with  local  time  stepping 
and  the  non-reflecting  boundary  condition  using  the  appropriate  0.  and  cui  for  global 
time  stepping  schemes.  The  result  is  that  the  convergence  behaviour  in  this  case  is 
unexpected.  In  fact,  the  convergence  is  better  than  with  local  time  stepping  with  the 
non-reflecting  condition  along  with  it’s  appropriate  and  wi.  The  bottom  dotted  line 
in  (9)  shows  the  convergence  behaviour  in  this  case,  and  another  order  of  magnitude 
is  obtained.  Figure  (10)  shows  the  pointwise  residuals  after  5000  iterations  in  this 
additional  case.  Clearly  the  best  steady  state  is  achieved  in  this  way. 

This  in  no  way  disproves  our  claim.  We  obtain  better  convergence  using  our 
non-reflecting  boundary  condition  than  with  extrapolation  conditions  both  in  global 
and  local  time  stepping  regimes.  However,  it  appears  that  our  method  of  finding  the 
coefficients  is  not  optimal,  so  that  perhaps  factors  other  than  long  wave  disturbances 
need  to  be  taken  into  consideration.  In  any  case,  the  important  result  that  we  have 
demonstrated  here  is  as  follows: 

Local  time  stepping  schemes  converge  to  steady  state  faster  than  global 
time  stepping  schemes  do  because  they  cause  longwave  disturbances  from 
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Figure  10;  Local  time  stepping  run  with  0  and  tui  taking  on  global  time  stepping 
values.  Residuals 


steady  state  to  decay  faster  than  schemes  utilizing  global  time  stepping. 


5.3  Compressible  Viscous  Flow  around  airfoils 

5.3.1  Case  1-  Global  time  stepping  without  multigrid  acceleration 

The  far-held  non-reflecting  boundary  condition  was  implemented  in  the  program  [9], 
as  a  means  to  test  its  effectiveness  in  the  compressible  viscous  subsonic  case  with 
Ri  =  5000  and  Moo  =  0.5.  An  appropriate  257  x  65  C-mesh  was  constructed  using 
a  hyperbolic  grid  generator.  Often,  hyperbolic  generators  create  grids  that  tend  to 
magnify  the  intrinsic  singularity  at  the  trailing  edge.  By  adding  additional  control 
points  we  have  overcome  this  phenomenon  and  generated  a  smooth  grid  reminiscent 
of  an  ellipticly  generated  one.  The  outflow  boundary  was  placed  5  chord  lengths 
away  from  the  trailing  edge,  and  was  made  as  perpendicular  as  possible  to  the  airfoil 
chord.  Two  runs  were  made,  one  with  extrapolation  boundary  conditions  and  the 
other  with  the  nonreflecting  boundary  condition.  In  each  run,  30000  time  steps  were 
taken.  Neither  multigrid  nor  other  accelerators  were  used.  Minimal  artificial  viscosity 
was  used  to  assure  convergence.  The  results  that  have  been  obtained  are  delineated 
below. 

The  and  norms  of  the  p  residual  are  shown  in  Figures  11  and  12-  Clearly 
the  nonreflecting  condition  accelerates  convergence  to  steady  state  more  than  the 
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Residuals  —  Maximum  Norm 


extrapolation  condition  does.  The  peaks  in  the  residuals  while  using  extrapolation 
have  not  yet  been  explained.  Due  to  stability  considerations  (long  time  numerical 
integration),  in  both  cases  the  residuals  increase  and  level  off  towards  the  end  of  the 
runs,  but  in  the  nonreflecting  case,  the  smallest  residuals  are  smaller  than  in  the 
extrapolation  case,  and  the  residuals  grow  sooner  and  faster  with  extrapolation  than 
with  the  nonreflecting  condition.  Hence,  the  nonreflecting  condition  produces  a  more 
robust  computational  environment,  and  achieves  a  better  steady  state. 

In  Figures  13  and  14,  the  pointwise  residuals  are  shown  for  extrapolation  and 
nonreflecting  boundary  conditions  for  30000  time  steps.  While  30000  time  steps  is 
clearly  beyond  the  optimal  computation  time,  it  is  still  evident  that  the  nonreflecting 
boundary  condition  is  more  effective  at  reducing  residuals  than  the  extrapolation 
condition.  In  addition,  in  Figure  15  the  values  of  p  and  u  along  the  grid  line:  j  =  5, 
i  =  I  —  258  are  shown.  This  grid  line  starts  in  the  wake  region  (far  field),  then  is  in 
the  boundary  layer  (around  the  airfoil),  and  returns  to  the  far  field.  Both  conditions 
give  identical  results  around  the  airfoil  -  in  the  boundary  layer.  However,  in  the  far 
field,  the  nonreflecting  condition  produces  a  smoother  solution  than  extrapolating 
does.  While  we  have  shown  this  result  for  j  =  5  only,  it  is  true  in  the  entire  wake 
region. 
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Average  Residuals 


Figure  12:  The  norm  of  the  p  residual 


5.3.2  Case  2-  Global  time  stepping  with  multigrid  acceleration 


In  order  to  decrease  the  residual  further  and  to  stabilize  the  runs,  we  have  investigated 
the  nonreflecting  boundary  condition  in  regimes  with  global  time  steps,  accelerated 
to  steady  state  with  multigrid  techniques. 

As  we  previously  mentioned,  the  nonreflecting  boundary  condition  needs  to  be 
modified  when  using  multigrid  acceleration  in  the  numerical  scheme  and  it  takes  the 
form  (4.6). 


In  this  case  we  compared  two  runs  -  one  run  using  multigrid  acceleration  and  ex¬ 
trapolation  outflow  conditions,  and  the  other  run  using  multigrid  acceleration  and  the 
above  nonreflecting  boundary  condition.  In  order  to  help  eliminate  short  wave  distur¬ 


bances  in  the  y  direction,  we  further  introduced  the  Fourier  smoothing  as  suggested 
by  Ryaben’kii  [8] 

1  .  1  „  1  „ 


+  +  iPl 


locally  near  the  lower  and  upper  corners  along  the  outflow  boundary  This  smoothing 
transformation  was  applied  each  time  our  nonreflecting  boundary  condition  was  called 
and  can  be  shown  to  identically  kill  off  the  shortest  wave  in  the  Fourier  expansion. 


In  the  figures  which  follow  (figures  16,  17  and  18),  we  show  the  convergence  be¬ 
havior  in  both  the  L°°  and  norms  in  using  at  outflow  extrapolation  and  the  nonre¬ 
flecting  condition  and  show  residual  plots  for  p  for  both  the  extrapolation  boundary 
condition  and  the  nonreflecting  boundary  condition. 
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The  same  C-grid  and  other  flow  parameters  were  used  in  these  calculations  as 
in  the  previous  subsection.  Multigrid  acceleration  was  applied  in  the  following  way: 
100  ly  cycles  (3  grids)  were  run  first  starting  on  the  first  coarser  grid  in  order  to 
set  up  a  reasonable  initial  condition.  The  results  of  these  sweeps  were  used  as  the 
initial  condition  for  a  set  of  2000  W  cycles  (4  grids)  starting  on  the  finest  grid.  Both 
second  order  and  fourth  order  artificial  viscosity  (see  [9]  and  [11])  were  applied  as 
well  as  implicit  residual  smoothing.  Residual  smoothing  insures  that  the  scheme  will 
converge  with  a  much  larger  CFL  number  —  hence  accelerating  convergence. 

The  results  are  as  follows:  The  residual  norms  obtained  using  both  boundary 
conditions  are  nearly  identical.  However  the  residual  plots  indicate  that  the  largest 
residuals  at  steady  state  with  the  both  outflow  conditions  are  at  the  corners  of  the 
outflow  boundary  (the  corner  problem)  and  at  inflow,  presumably  propagated  back¬ 
ward  from  the  corners.  The  maximum  norm,  L°°  of  the  steady  state  residual  while 
using  the  nonreflecting  boundary  condition  at  outflow,  is  slightly  greater  then  when 
using  the  extrapolation  conditions  at  outflow  as  is  confirmed  in  figure  16.  However, 
in  the  wake  region,  the  residuals  are  smaller  when  using  the  nonreflecting  boundary 
condition  then  when  using  extrapolation  conditions.  This  is  clearly  seen  in  figures  17 
and  18. 
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Figure  14:  Residuals  using  Nonreflecting  conditions  for  (a)  p,  (b)  pu,  (c)  pv  and 
(d)  e 

6  Conclusion 


We  have  described  the  derivation  of  a  nonreflecting  far  field  boundary  condition  for 
steady  viscous,  compressible,  two  dimensional  external  flows.  This  condition  pro¬ 
duces  a  smoother  far  field  solution  than  the  presently  used  extrapolation  boundary 
conditions,  while  retaining  the  solution  in  the  boundary  layer.  It  also  smooths  out 
the  trailing  edge  singularity  by  one  to  two  orders  of  magnitude.  The  method  can 
be  applied  to  incompressible  flows,  where  it  is  just  as  effective  in  driving  residuals 
down  to  steady  state  faster  than  extrapolation  conditions.  In  addition,  nonreflecting 
boundary  conditions  of  the  type  described  in  this  work  can  be  developed  for  schemes 
whose  time  stepping  is  local,  as  well  as  for  schemes  utilizing  multigrid  acceleration. 

In  addition  to  the  benefits  accrued  from  using  out  boundary  condition,  it  is  easy 
to  program,  and  does  not  require  the  use  of  much  additional  memory  in  an  existing 
program.  Also,  its  simplicity  and  local  form  require  little  cpu  time. 

Perhaps  the  most  important  feature  of  our  approach  is  that  while  in  this  paper  we 
have  chosen  to  deal  with  steady  state  problems,  our  methodology  could  be  used  to 
develop  a  condition  for  evolution  equations  that  do  not  reach  a  steady  state.  In  this 
case,  we  would  perturb  the  linerized  equations  about  the  appropriate  wave  group  that 
we  wish  to  absorb  (not  necessarily  about  long  waves)  and  proceed  exactly  as  we  did 
for  steady  state  flows.  Of  course  an  appropriate  far  field  profile  needs  to  be  formulated 
as  well  for  each  flow,  which  now  may  be  time  dependent.  Hence  our  methodolgy  is 
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Figure  15:  The  values  of  p  and  u  along  the  grid  line:  j  =  5,  i  =  1  —  258 


versatile  and  is  shown  to  be  effective  in  reducing  residuals  in  external  viscous  flows. 
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Figure  17:  The  p  residuals  after  2000  multigrid  W  cycles.  Extrapolation  outflow 
conditions. 


Density  Residuals-  Multigrid 


Figure  18:  The  p  residuals  after  2000  multigrid  W  cycles.  Nonreflecting  outflow 
conditions. 
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